Suplementary material: Generalized method of moments and identifiability

knitr::opts_chunk$set(echo = TRUE)
library(latex2exp)
library(IsoPriceR)

Attaching package: ‘IsoPriceR’

The following object is masked _by_ ‘.GlobalEnv’:

    sample_g

library(ggplot2)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
getwd()

[1] "/Users/goffard/code/market_based_insurance_ratemaking/notebooks"

Let \(X\) be a random variable that represents the total health related over a given time period, say one year, associated to a pet. Our goal is to find the parameter \(\theta\) that best explains our market data made of insurance quotes \[ \tilde{p}_i = f_i\left[\mathbb{E}(g_i(X))\right],\text{ }i = 1,\ldots, n, \] where \(g_i\) are the coverage functions and \(f_i\) are the loading function. The coverage functions are known and of the form \[ g(x) = \min(\max(r\cdot x - d, 0), l), \] where \(r\) is the rate of coverage, \(d\) is the deductible and \(l\) is the limit.The loading functions are unknown and will be approximated by a generic function \(f\).

In the first section we consider the problem of finding \(\theta\) if we know access to the pure premiums

\[ p_i = \left[\mathbb{E}(g_i(X))\right],\text{ }i = 1,\ldots, n, \] which is not a realistic situation in practice. We then consider an isotonic regression model to approximate \(f\) and compares it to a parametric curve, linear withe \(f(x) = a + b\cdot x\) for simplicity. We consider an example where

\[ X = \sum_{k=1}^NU_k, \] where \(N\sim\text{Pois}(\lambda = 3)\) and the \(U_k\) are iid and lognormally distributed \(U\sim\text{LogNormal}(\mu =0 , \sigma = 1)\)

set.seed(142)

# Loss model
freq_dist <- model("poisson", c("lambda"))
sev_dist <- model("lognormal", c("mu", "sigma"))
l_m <- loss_model(freq_dist, sev_dist, "Poisson-Lognormal")

# True parameter
lambda_true <- 3; mu_true <- 0; sigma_true <- 1
true_parms <- t(as.matrix(c(lambda_true, mu_true, sigma_true)))
colnames(true_parms) <- l_m$parm_names

# Data
## Sample of claim data to calculate the true pure premium via Monte Carlo
X_true <- sample_X(l_m, true_parms, R = 10000)
summary(X_true)

       V1        
 Min.   : 0.000  
 1st Qu.: 1.705  
 Median : 3.765  
 Mean   : 4.946  
 3rd Qu.: 6.772  
 Max.   :76.830  

I - Generalized method of moments

If the pure premiums

\[ p_i = \mathbb{E}_(g_i(X)),\text{ }i = 1,\ldots, n, \]

are available then our estimation problem becomes find \(\theta\) that minimizes

\[ d\left[p_{1:n},p_{1:n}^\theta\right] = \sum_i w_i^{RMSE}(p_i -p_i^\theta)^2, \]

where \(p_i^\theta = \mathbb{E}_{\theta}(g_i(X))\). Identifiability holds if there exists a unique \(\theta^*\) such that

\[ \theta^* = \underset{\theta\in \Theta}{\text{argmin}}\, d\left[p_{1:n},p_{1:n}^\theta\right]. \]

This condition ressembles the global identification conditions of the Generalized Method of Moments which is typically hard to show. Existence of minimum requires that \(\theta\mapsto d\left[p_{1:n},p_{1:n}^\theta\right]\) is continuous and that \(\Theta\) is compact. Uniqueness would then hold if \(\theta\mapsto d\left[p_{1:n},p_{1:n}^\theta\right]\) is strictly convex. This condition is difficult to verify because an analytical expression of \(\theta\mapsto d\left[p_{1:n},p_{1:n}^\theta\right]\) is out of our reach and this for almost all the compound sum \(X\). A necessary condition is that the dimension of \(\theta\) is lower than \(n\) which is the number of moments we hold. Let us consider the case \(n = 3\) with

\[ g_1 = 0.85\cdot X\text{, } g_2 = \max(X - 1.8, 0) \text{, and }g_3 = \min(X, 6) \]

it is equivalent to the coverages \((r_1 = 0.85, d_1 = 0, l_1 = \infty)\), \((r_2 = 1, d_2 = 1.8, l_2 = \infty)\), and \((r_3 = 1, d_3 = 0, l_3 = 6)\).

## Characteristics of the pure premium which depends on the parameter of the insurance coverage
p_p <- pure_premium("xs", c("r", "d", "l"))
## Random premium parameters
n <- 3; r <- c(0.75, 1, 1); d <- c(0, 1.8, 0); l <- c(Inf, Inf, 6)

## Data frame with the insurance coverage parameters and the pure premium
th_premium <- matrix(c(r, d, l), nrow = n)
pps <- compute_pure_premia(X_true, coverage_type = "xs", th_premium)
sds <- sd_pure_premia(X_true, coverage_type = "xs", th_premium)
df_pp <- data.frame(r = r, d = d, l = l, pp=pps)
df_pp

The pure premium computed above corresponds to a well specified case where there exists \(\theta^* = \theta_0\) for which

\[ d\left[p_{1:n},p_{1:n}^{\theta^{\ast}}\right] = 0 \] Note that we assign the same weight to each if the pre premium. We can check the uniqueness empirically, first in one dimension as we plot the functions \(\sigma\mapsto d\left[p_{1:n},p_{1:n}^{\theta^{\ast}}\right]\Big\rvert_{\mu =0, \lambda = 3}\), \(\mu\mapsto d\left[p_{1:n},p_{1:n}^{\theta^{\ast}}\right]\Big\rvert_{\sigma =1, \lambda = 3}\), and \(\lambda\mapsto d\left[p_{1:n},p_{1:n}^{\theta^{\ast}}\right]\Big\rvert_{\mu =0, \sigma = 1}\)

data <- df_pp[,c("r", "d", "l")]
data["x"] <- df_pp$pp
ths_premium <- as.matrix(data[p_p$parm_names])
lambda_true <- 3; mu_true <- 0; sigma_true <- 1

# Sigma --------------
# Pure premium calculation for a range of risk model parameters
sigmas <- seq(0.1, 2, 0.1)
cloud <- as.matrix(expand.grid(c(3), c(0),sigmas))
colnames(cloud) <- l_m$parm_names
X <- sample_X(l_m , cloud, R = 10000)
pp_fake <- compute_pure_premia(X, coverage_type = "xs", ths_premium)
                                                                     
# Calculation of the RMSE between true and wrong pure premium
ds <- sapply(1:nrow(cloud), function(k) mean((pp_fake[, k] - data$x)^2)^(1/2)) 
df_ds <- data.frame(cloud)
df_ds <- cbind(df_ds, ds)
(sig_pure_RMSE <- ggplot(data = df_ds) + geom_line(mapping = aes(x = sigma, y = ds)) +  ggplot2::theme_classic(base_size = 14) + labs(x = expression(sigma), title = "RMSE", y = "") )
ggsave("../figures/sig_pure_RMSE.pdf", sig_pure_RMSE)

Saving 7.29 x 4.51 in image

# Lambda --------------
lambdas <- seq(0.1, 5, 0.1)
cloud <- as.matrix(expand.grid(lambdas, c(0),c(1)))
colnames(cloud) <- l_m$parm_names
X <- sample_X(l_m , cloud, R = 10000)
pp_fake <- compute_pure_premia(X, coverage_type = "xs", ths_premium)
                                                                     
# Calculation of the RMSE between true and wrong pure premium
ds <- sapply(1:nrow(cloud), function(k) mean((pp_fake[, k] - data$x)^2)^(1/2)) 
df_ds <- data.frame(cloud)
df_ds <- cbind(df_ds, ds)
(lambda_pure_RMSE <- ggplot(data = df_ds) + geom_line(mapping = aes(x = lambda, y = ds)) +  ggplot2::theme_classic(base_size = 14) + labs(x = expression(lambda), title = "RMSE", y = "") )
ggsave("../figures/lambda_pure_RMSE.pdf", lambda_pure_RMSE)

Saving 7.29 x 4.51 in image

# mu --------------
mus <- seq(-3, 3, 0.1)
cloud <- as.matrix(expand.grid(c(3), mus,c(1)))
colnames(cloud) <- l_m$parm_names
X <- sample_X(l_m , cloud, R = 10000)
pp_fake <- compute_pure_premia(X, coverage_type = "xs", ths_premium)
                                                                     
# Calculation of the RMSE between true and wrong pure premium
ds <- sapply(1:nrow(cloud), function(k) mean((pp_fake[, k] - data$x)^2)^(1/2)) 
df_ds <- data.frame(cloud)
df_ds <- cbind(df_ds, ds)
(mu_pure_RMSE <- ggplot(data = df_ds) + geom_line(mapping = aes(x = mu, y = ds)) +  ggplot2::theme_classic(base_size = 14) + labs(x = expression(mu), title = "RMSE", y = "") )
ggsave("../figures/mu_pure_RMSE.pdf", mu_pure_RMSE)

Saving 7.29 x 4.51 in image

Each parameters can be identified separately, now of course identifying them jointly is another matter, there might exists several combinations of parameters that yield the same pure premium. Our optimization allows us to investigate by searching efficiently the parameter space.

We use the follwing prior assumptions: \[ \lambda\sim\text{Unif}(0, 10)\text{, }\mu\sim\text{Unif}(-3, 3)\text{, and }\sigma\sim\text{Unif}(0, 5) \] and the following settings for the algorithm

\[ J = 500\text{ and }R = 500 \] Recall that \(J\) is the population size of the clouds of particles and \(R\) is the number of Monte Carlo replications.

# Prior settings for the parameters
prior_lambda <- prior_dist("uniform", "lambda",c(0, 10))
prior_mu <- prior_dist("uniform", "mu",c(-3,3))
prior_sigma <- prior_dist("uniform", "sigma", c(0, 5))
prior_distributions <- list(prior_lambda, prior_mu, prior_sigma)
model_prior <- independent_priors(prior_distributions)

# Settings for the abc algorithm
popSize <- 500;  MC_R<- 500; acc <- 10; sp_bounds <- c(1, 1); eps_min <- 0.02

sum(sds/sqrt(MC_R)^2)

data <- df_pp[,c("r", "d", "l")]
data["x"] <- df_pp$pp
# res_abc <- abc(data, model_prior, l_m, p_p, popSize, MC_R, acc, sp_bounds, eps_min)
# save(res_abc, file = "../data/res_abc_supp_material_Section1.RData")

load("../data/res_abc_supp_material_Section1.RData")
# We retrieve the cloud of particles associated to each generation
cloud_dfs <- lapply(res_abc$Clouds, as.data.frame)
cloud_df <- do.call(rbind, cloud_dfs)
cloud_df$generation <- sort(rep(1:length(res_abc$Clouds), popSize))
cloud_df_plot <- rbind(
  data.frame(parm = cloud_df$lambda, parm_name = "lambda", generation = cloud_df$generation),
  data.frame(parm = cloud_df$sigma, parm_name = "sigma", generation = cloud_df$generation)
)
cloud_df_plot$parm_name <-as.factor(cloud_df_plot$parm_name) 
levels(cloud_df_plot$parm_name) <- c(lambda = TeX("$\\lambda$"), sigma = TeX("$\\sigma$"))

(posterior_lambda <- ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=lambda, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\lambda$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,1]), color="black",
                                   linetype="dashed"))
ggsave("../figures/posterior_lambda_supp_mat.pdf", posterior_lambda)

(posterior_mu <- ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=mu, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\mu$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,2]), color="black",
                                   linetype="dashed"))
ggsave("../figures/posterior_mu_supp_mat.pdf", posterior_mu)

(posterior_sig <- ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=sigma, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\sigma$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,3]), color="black",
                                   linetype="dashed"))
ggsave("../figures/posterior_sig_supp_mat.pdf", posterior_sig)

II - Optimization based on the commercial premiums without regularization

We do not have access to the pure premium, instead we have the commercial premiums

\[ \tilde{p}_i = f_i\left[\mathbb{E}(g_i(X))\right],\text{ }i = 1,\ldots, n, \]

where the loading functions must be approximated. We therefore have to choose a model \(f\) such that

\[ \tilde{p}_i = f(p_{i}^\theta),\text{ }i = 1,\ldots, n \] and hope that there exists a unique \(\theta^*\) such that

\[ \theta^* = \underset{\theta\in \Theta}{\text{argmin}}\, d\left[\tilde{p}_{1:n},f(p_{1:n}^\theta)\right]. \] This typically depends on the choice for \(f\) we consider isotonic and linear. Isotonic regression impose monotonicity which is desirable since we believe it would make sense if

\[ p_i > p_j\Rightarrow \tilde{p}_i > \tilde{p}_j. \]

Furthermore isotonic regression is non parametric which mitigate the risk of misspecification. The drawback is overfitting as perfect interpolation of may occur in a small sample size setting.

Well-specified case

In our example we take the following loading fuunctions

\[ f_1(x) = 1.1\cdot x\text{, }f_2(x) = 1.2\cdot x \text{, and } f_3(x) = 1.3\cdot x \]

data <- df_pp[,c("r", "d", "l")]
data["x"] <- c(df_pp$pp[1]*1.38 , df_pp$pp[2]*1.1 , df_pp$pp[3] * 1.4 )
# data["x"] <- df_pp$pp
ths_premium <- as.matrix(data[p_p$parm_names])
lambda_true <- 3; mu_true <- 0; sigma_true <- 1
(pp_cp_well_specified <- ggplot() + geom_point(mapping = aes(x = df_pp$pp, y = data$x)) +  ggplot2::theme_classic(base_size = 14) + labs(x = "Pure Premium", title = "Commercial premium", y = "") )
ggsave("../figures/pp_cp_well_specified.pdf", pp_cp_well_specified)

Saving 7.29 x 4.51 in image

# Sigma --------------
# Pure premium calculation for a range of risk model parameters
sigmas <- seq(0.1, 2, 0.02)
cloud <- as.matrix(expand.grid(c(3), c(0),sigmas))
colnames(cloud) <- l_m$parm_names
X <- sample_X(l_m , cloud, R = 10000)
pp_fake <- compute_pure_premia(X, coverage_type = "xs", ths_premium)
# Calculation of the RMSE between true and wrong commercial premium
ress_iso <- lapply(1:ncol(pp_fake), function(k) isoreg(pp_fake[, k], data$x))

ress_lm <- sapply(1:ncol(pp_fake), function(k) as.numeric(lm(data$x~pp_fake[, k] - 1)$coefficients))
ress_lm[sigmas == 1]

[1] 1.297967

as.numeric(lm(data$x ~ pp_fake[, sigmas == 1] - 1)$coefficients)

[1] 1.297967

mean((mean((pp_fake[, sigmas == 1] * ress_lm[sigmas == 1] - data$x)^2)^(1/2) - data$x)^2)^(1/2)

[1] 4.20526

fs_iso <- lapply(ress_iso, function(res_iso) stepfun(sort(res_iso$x), c(min(res_iso$yf) /2 ,
                                                                      res_iso$yf), f = 1))
ds_iso <- sapply(1:nrow(cloud), function(k) mean((fs_iso[[k]](pp_fake[, k]) - data$x)^2)^(1/2)) 
ds_lm <- sapply(1:nrow(cloud), function(k) mean((pp_fake[, k] / ress_lm[k] - data$x)^2)^(1/2)) 

df_ds <- data.frame(cloud)
df_ds <- cbind(df_ds, ds_iso)
df_ds <- cbind(df_ds, ds_lm)
ggplot(data = df_ds) + geom_line(mapping = aes(x = sigma, y = ds_iso))

ggplot(data = df_ds) + geom_line(mapping = aes(x = sigma, y = ds_lm))

sigmas[ds_lm == min(ds_lm)]

[1] 1.22

# Prior settings for the parameters
prior_lambda <- prior_dist("uniform", "lambda",c(0, 5))
prior_mu <- prior_dist("uniform", "mu",c(-3,3))
prior_sigma <- prior_dist("uniform", "sigma", c(0, 2))
prior_distributions <- list(prior_lambda, prior_mu, prior_sigma)
model_prior <- independent_priors(prior_distributions)

# Settings for the abc algorithm
popSize <- 500;  MC_R<- 500; acc <- 0.01; sp_bounds <- c(10^(-16), Inf); eps_min <- 0.01
res_abc <- abc(data, model_prior, l_m, p_p, popSize, MC_R, acc, sp_bounds, eps_min)

[1] "Sampling generation of particles  1 ; eps =  Inf ; ESS =  500"
[1] "Sampling generation of particles  2 ; eps =  0.638654634426981 ; ESS =  250"
===============================================================================================================================================================================================
[1] 23
[1] "Sampling generation of particles  3 ; eps =  0.0179453779610254 ; ESS =  250.653519160833"
===============================================================================================================================================================================================
[1] 500
[1] "Sampling generation of particles  4 ; eps =  0.01 ; ESS =  250.706714954091"
===============================================================================================================================================================================================
[1] 500

# We retrieve the cloud of particles associated to each generation
cloud_dfs <- lapply(res_abc$Clouds, as.data.frame)
cloud_df <- do.call(rbind, cloud_dfs)
cloud_df$generation <- sort(rep(1:length(res_abc$Clouds), popSize))
cloud_df_plot <- rbind(
  data.frame(parm = cloud_df$lambda, parm_name = "lambda", generation = cloud_df$generation),
  data.frame(parm = cloud_df$sigma, parm_name = "sigma", generation = cloud_df$generation)
)
cloud_df_plot$parm_name <-as.factor(cloud_df_plot$parm_name) 
levels(cloud_df_plot$parm_name) <- c(lambda = TeX("$\\lambda$"), sigma = TeX("$\\sigma$"))
ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=lambda, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\lambda$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,1]), color="black",
                                   linetype="dashed")

ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=mu, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\mu$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,2]), color="black",
                                   linetype="dashed")

ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=sigma, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\sigma$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,3]), color="black",
                                   linetype="dashed")

III - Optimization based on the commercial premiums with regularization

# Adding the regularization terms
sp_bounds <- c(min(df_pp$pp / data$x), max(df_pp$pp / data$x))
# sp_bounds <- c(0.6, 0.95)
ds1 <- sapply(1:ncol(pp_fake), function(k) sqrt(mean((pmax(data$x - pp_fake[,k] /sp_bounds[1], 0) )^2)))
ds2 <- sapply(1:ncol(pp_fake), function(k) sqrt(mean((pmax( pp_fake[,k] / sp_bounds[2] - data$x , 0) )^2)))
ds_iso_reg <- ds_iso + ds1 + ds2
ds_lm_reg <- ds_lm + ds1 + ds2
df_ds <- data.frame(cloud)
df_ds <- cbind(df_ds, ds_iso_reg, ds_lm_reg)
ggplot(data = df_ds) + geom_line(mapping = aes(x = sigma, y = ds_iso_reg))

ggplot(data = df_ds) + geom_line(mapping = aes(x = sigma, y = ds_lm_reg))

print(sigmas[ds_lm_reg==min(ds_lm_reg)])

[1] 1.16

print(sigmas[ds_iso_reg==min(ds_iso_reg)])

[1] 1

# Prior settings for the parameters
prior_lambda <- prior_dist("uniform", "lambda",c(0, 5))
prior_mu <- prior_dist("uniform", "mu",c(-3,3))
prior_sigma <- prior_dist("uniform", "sigma", c(0, 2))
prior_distributions <- list(prior_lambda, prior_mu, prior_sigma)
model_prior <- independent_priors(prior_distributions)

# Settings for the abc algorithm
popSize <- 500;  MC_R<- 500; acc <- 0.001; sp_bounds <- c(min(df_pp$pp / data$x), max(df_pp$pp / data$x))
eps_min <- 0.0001
res_abc <- abc(data, model_prior, l_m, p_p, popSize, MC_R, acc, sp_bounds, eps_min)

[1] "Sampling generation of particles  1 ; eps =  Inf ; ESS =  500"
[1] "Sampling generation of particles  2 ; eps =  4.43689039270973 ; ESS =  250"
===============================================================================================================================================================================================
[1] 260
[1] "Sampling generation of particles  3 ; eps =  3.35330585831361 ; ESS =  250.772464133936"
===============================================================================================================================================================================================
[1] 263
[1] "Sampling generation of particles  4 ; eps =  2.23229920434542 ; ESS =  250.886475662332"
===============================================================================================================================================================================================
[1] 257
[1] "Sampling generation of particles  5 ; eps =  1.34653967981587 ; ESS =  250.030381611654"
===============================================================================================================================================================================================
[1] 259
[1] "Sampling generation of particles  6 ; eps =  0.731400077523843 ; ESS =  250.89930191337"
===============================================================================================================================================================================================
[1] 258
[1] "Sampling generation of particles  7 ; eps =  0.319613497124356 ; ESS =  250.2391400301"
===============================================================================================================================================================================================
[1] 260
[1] "Sampling generation of particles  8 ; eps =  0.148285269745679 ; ESS =  250.259059634864"
===============================================================================================================================================================================================
[1] 268
[1] "Sampling generation of particles  9 ; eps =  0.0739579668388663 ; ESS =  250.0545936104"
===============================================================================================================================================================================================
[1] 268
[1] "Sampling generation of particles  10 ; eps =  0.021013252054564 ; ESS =  250.599302376695"
===============================================================================================================================================================================================
[1] 15
[1] "Sampling generation of particles  11 ; eps =  0.0179453779610254 ; ESS =  250.574578198601"
===============================================================================================================================================================================================
[1] 290
[1] "Sampling generation of particles  12 ; eps =  0.00378496163232778 ; ESS =  250.632080483012"
===============================================================================================================================================================================================
[1] 427
[1] "Sampling generation of particles  13 ; eps =  1e-04 ; ESS =  250.70095264432"
===============================================================================================================================================================================================
[1] 500
[1] "Sampling generation of particles  14 ; eps =  1e-04 ; ESS =  250.701552548448"
===============================================================================================================================================================================================
[1] 500

# We retrieve the cloud of particles associated to each generation
cloud_dfs <- lapply(res_abc$Clouds, as.data.frame)
cloud_df <- do.call(rbind, cloud_dfs)
cloud_df$generation <- sort(rep(1:length(res_abc$Clouds), popSize))
cloud_df_plot <- rbind(
  data.frame(parm = cloud_df$lambda, parm_name = "lambda", generation = cloud_df$generation),
  data.frame(parm = cloud_df$sigma, parm_name = "sigma", generation = cloud_df$generation)
)
cloud_df_plot$parm_name <-as.factor(cloud_df_plot$parm_name) 
levels(cloud_df_plot$parm_name) <- c(lambda = TeX("$\\lambda$"), sigma = TeX("$\\sigma$"))
ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=lambda, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\lambda$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,1]), color="black",
                                   linetype="dashed")

ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=mu, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\mu$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,2]), color="black",
                                   linetype="dashed")

ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=sigma, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\sigma$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,3]), color="black",
                                   linetype="dashed")

res_abc$ds[[13]]

  [1] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
 [15] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
 [29] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 8.540344e-05 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
 [43] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
 [57] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
 [71] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
 [85] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
 [99] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[113] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[127] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[141] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[155] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[169] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[183] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[197] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[211] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[225] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[239] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[253] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[267] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[281] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[295] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[309] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[323] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[337] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[351] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[365] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[379] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[393] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[407] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[421] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[435] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[449] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[463] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[477] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16
[491] 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16 7.251946e-16

Well-specified case

In our example we take the following loading fuunctions

\[ f_1(x) = 1.1\cdot x\text{, }f_2(x) = 1.2\cdot x \text{, and } f_3(x) = 1.3\cdot x \]

data <- df_pp[,c("r", "d", "l")]
data["x"] <- c(df_pp$pp[1]*1.1 , df_pp$pp[2]*1.2 , df_pp$pp[3] * 1.3 )
# data["x"] <- df_pp$pp
ths_premium <- as.matrix(data[p_p$parm_names])
lambda_true <- 3; mu_true <- 0; sigma_true <- 1
(pp_cp_well_specified <- ggplot() + geom_point(mapping = aes(x = df_pp$pp, y = data$x)) +  ggplot2::theme_classic(base_size = 14) + labs(x = "Pure Premium", title = "Commercial premium", y = "") )
ggsave("../figures/pp_cp_well_specified.pdf", pp_cp_well_specified)

Saving 7.29 x 4.51 in image

# Sigma --------------
# Pure premium calculation for a range of risk model parameters
sigmas <- seq(0.1, 2, 0.02)
cloud <- as.matrix(expand.grid(c(3), c(0),sigmas))
colnames(cloud) <- l_m$parm_names
X <- sample_X(l_m , cloud, R = 10000)
pp_fake <- compute_pure_premia(X, coverage_type = "xs", ths_premium)
# Calculation of the RMSE between true and wrong commercial premium
ress_iso <- lapply(1:ncol(pp_fake), function(k) isoreg(pp_fake[, k], data$x))

ress_lm <- sapply(1:ncol(pp_fake), function(k) as.numeric(lm(data$x~pp_fake[, k] - 1)$coefficients))
ress_lm[sigmas == 1]

[1] 1.19087

as.numeric(lm(data$x ~ pp_fake[, sigmas == 1] - 1)$coefficients)

[1] 1.19087

mean((mean((pp_fake[, sigmas == 1] * ress_lm[sigmas == 1] - data$x)^2)^(1/2) - data$x)^2)^(1/2)

[1] 3.984624

fs_iso <- lapply(ress_iso, function(res_iso) stepfun(sort(res_iso$x), c(min(res_iso$yf) /2 ,
                                                                      res_iso$yf), f = 1))
ds_iso <- sapply(1:nrow(cloud), function(k) mean((fs_iso[[k]](pp_fake[, k]) - data$x)^2)^(1/2)) 
ds_lm <- sapply(1:nrow(cloud), function(k) mean((pp_fake[, k] / ress_lm[k] - data$x)^2)^(1/2)) 

df_ds <- data.frame(cloud)
df_ds <- cbind(df_ds, ds_iso)
df_ds <- cbind(df_ds, ds_lm)
ggplot(data = df_ds) + geom_line(mapping = aes(x = sigma, y = ds_iso))

ggplot(data = df_ds) + geom_line(mapping = aes(x = sigma, y = ds_lm))

sigmas[ds_lm == min(ds_lm)]

[1] 1.16

# Prior settings for the parameters
prior_lambda <- prior_dist("uniform", "lambda",c(0, 5))
prior_mu <- prior_dist("uniform", "mu",c(-3,3))
prior_sigma <- prior_dist("uniform", "sigma", c(0, 2))
prior_distributions <- list(prior_lambda, prior_mu, prior_sigma)
model_prior <- independent_priors(prior_distributions)

# Settings for the abc algorithm
popSize <- 500;  MC_R<- 500; acc <- 0.01; sp_bounds <- c(10^(-16), Inf); eps_min <- 0.01
res_abc <- abc(data, model_prior, l_m, p_p, popSize, MC_R, acc, sp_bounds, eps_min)

[1] "Sampling generation of particles  1 ; eps =  Inf ; ESS =  500"
[1] "Sampling generation of particles  2 ; eps =  0.297386732987592 ; ESS =  250"
===============================================================================================================================================================================================
[1] 242
[1] "Sampling generation of particles  3 ; eps =  0.258092879599686 ; ESS =  250.503338748765"
===============================================================================================================================================================================================
[1] 500
[1] "Sampling generation of particles  4 ; eps =  0.01 ; ESS =  250.445785786703"
===============================================================================================================================================================================================
[1] 500
[1] "Sampling generation of particles  5 ; eps =  0.01 ; ESS =  250.653120238212"
===============================================================================================================================================================================================
[1] 500

# We retrieve the cloud of particles associated to each generation
cloud_dfs <- lapply(res_abc$Clouds, as.data.frame)
cloud_df <- do.call(rbind, cloud_dfs)
cloud_df$generation <- sort(rep(1:length(res_abc$Clouds), popSize))
cloud_df_plot <- rbind(
  data.frame(parm = cloud_df$lambda, parm_name = "lambda", generation = cloud_df$generation),
  data.frame(parm = cloud_df$sigma, parm_name = "sigma", generation = cloud_df$generation)
)
cloud_df_plot$parm_name <-as.factor(cloud_df_plot$parm_name) 
levels(cloud_df_plot$parm_name) <- c(lambda = TeX("$\\lambda$"), sigma = TeX("$\\sigma$"))
ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=lambda, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\lambda$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,1]), color="black",
                                   linetype="dashed")

ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=mu, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\mu$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,2]), color="black",
                                   linetype="dashed")

ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=sigma, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\sigma$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,3]), color="black",
                                   linetype="dashed")

III - Optimization based on the commercial premiums with regularization

print(sigmas[ds_iso_reg==min(ds_iso_reg)]); print(min(ds_iso_reg))

[1] 0.96
[1] 0.001560745

# Prior settings for the parameters
prior_lambda <- prior_dist("uniform", "lambda",c(0, 5))
prior_mu <- prior_dist("uniform", "mu",c(-3,3))
prior_sigma <- prior_dist("uniform", "sigma", c(0, 2))
prior_distributions <- list(prior_lambda, prior_mu, prior_sigma)
model_prior <- independent_priors(prior_distributions)

# Settings for the abc algorithm
popSize <- 500;  MC_R<- 500; acc <- 0.002; sp_bounds <- c(min(df_pp$pp / data$x), max(df_pp$pp / data$x))
eps_min <- 0.002
res_abc <- abc(data, model_prior, l_m, p_p, popSize, MC_R, acc, sp_bounds, eps_min)

[1] "Sampling generation of particles  1 ; eps =  Inf ; ESS =  500"
[1] "Sampling generation of particles  2 ; eps =  3.94292337790187 ; ESS =  250"
===============================================================================================================================================================================================
[1] 264
[1] "Sampling generation of particles  3 ; eps =  2.88561559648946 ; ESS =  250.910601767801"
===============================================================================================================================================================================================
[1] 261
[1] "Sampling generation of particles  4 ; eps =  1.82835959111493 ; ESS =  250.88826264432"
===============================================================================================================================================================================================
[1] 257
[1] "Sampling generation of particles  5 ; eps =  1.12298655265911 ; ESS =  250.516480072063"
===============================================================================================================================================================================================
[1] 260
[1] "Sampling generation of particles  6 ; eps =  0.674098340085706 ; ESS =  250.233588685872"
===============================================================================================================================================================================================
[1] 256
[1] "Sampling generation of particles  7 ; eps =  0.350334206408645 ; ESS =  250.129029112909"
===============================================================================================================================================================================================
[1] 259
[1] "Sampling generation of particles  8 ; eps =  0.134850428265648 ; ESS =  250.06467485184"
===============================================================================================================================================================================================
[1] 267
[1] "Sampling generation of particles  9 ; eps =  0.0537250475762552 ; ESS =  250.750732391304"
===============================================================================================================================================================================================
[1] 270
[1] "Sampling generation of particles  10 ; eps =  0.00907367331436939 ; ESS =  250.332252647683"
===============================================================================================================================================================================================
[1] 438
[1] "Sampling generation of particles  11 ; eps =  0.002 ; ESS =  250.655511831248"
===============================================================================================================================================================================================
[1] 500
[1] "Sampling generation of particles  12 ; eps =  0.002 ; ESS =  250.455240112145"
===============================================================================================================================================================================================
[1] 500

# We retrieve the cloud of particles associated to each generation
cloud_dfs <- lapply(res_abc$Clouds, as.data.frame)
cloud_df <- do.call(rbind, cloud_dfs)
cloud_df$generation <- sort(rep(1:length(res_abc$Clouds), popSize))
cloud_df_plot <- rbind(
  data.frame(parm = cloud_df$lambda, parm_name = "lambda", generation = cloud_df$generation),
  data.frame(parm = cloud_df$sigma, parm_name = "sigma", generation = cloud_df$generation)
)
cloud_df_plot$parm_name <-as.factor(cloud_df_plot$parm_name) 
levels(cloud_df_plot$parm_name) <- c(lambda = TeX("$\\lambda$"), sigma = TeX("$\\sigma$"))
ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=lambda, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\lambda$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,1]), color="black",
                                   linetype="dashed")

ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=mu, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\mu$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,2]), color="black",
                                   linetype="dashed")

ggplot2::ggplot(data = cloud_df ) + ggplot2::geom_density(data = cloud_df, mapping = aes(x=sigma, color = as.factor(generation))) + ggplot2::labs(x=TeX("$\\sigma$"), y = "Density", color = "Generation") + ggplot2::theme_classic(base_size = 20) + ggplot2::theme(legend.position="right") + ggplot2::geom_vline(ggplot2::aes(xintercept=true_parms[,3]), color="black",
                                   linetype="dashed")

res_abc$ds[[12]]

  [1] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
 [16] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
 [31] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
 [46] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
 [61] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
 [76] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
 [91] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[106] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[121] 0.001100374 0.001100374 0.001100374 0.001102460 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[136] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[151] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[166] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001842784 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[181] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[196] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[211] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[226] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[241] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001900111
[256] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001192304 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[271] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[286] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[301] 0.001100374 0.001100374 0.001108272 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[316] 0.001100374 0.001155256 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[331] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[346] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[361] 0.001100374 0.001100374 0.001100374 0.001422143 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[376] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001255541 0.001100374
[391] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001338872 0.001100374 0.001100374 0.001100374 0.001100374
[406] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[421] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[436] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[451] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[466] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[481] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374
[496] 0.001100374 0.001100374 0.001100374 0.001100374 0.001100374